﻿using MathNet.Numerics.LinearAlgebra.Factorization;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearRegression;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace UtilZ.Dotnet.MathNetEx.LinearAlgebra
{
    /// <summary>
    /// 
    /// </summary>
    public class FitEx
    {
        /// <summary>
        /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k,
        /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array, compatible with Polynomial.Evaluate.
        /// A polynomial with order/degree k has (k+1) coefficients and thus requires at least (k+1) samples.
        /// </summary>
        public static double[] Polynomial(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR)
        {
            //MatrixBuilder<double> builder = Matrix<double>.Build;
            //MatrixBuilder<double> builder = BuilderInstance<double>.Matrix;
            //var builder = new MathNet.Numerics.LinearAlgebra.Double.MatrixBuilder();
            var storage = MathNet.Numerics.LinearAlgebra.Storage.DenseColumnMajorMatrixStorage<double>.OfInit(x.Length, order + 1, (i, j) => Math.Pow(x[i], j));
            var design = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(storage);
            var designCopy = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(storage);
            //var design = builder.Dense(x.Length, order + 1, (i, j) => Math.Pow(x[i], j));

            VectorBuilder<double> builder2 = DoubleVectorBuilderEx.Builder;
            Vector<double> yy = builder2.Dense(y);
            //return MultipleRegression.DirectMethod(design, yy, method).ToArray();

            if (design.RowCount != yy.Count)
            {
                throw new ArgumentException($"All sample vectors must have the same length. However, vectors with disagreeing length {design.RowCount} and {yy.Count} have been provided. A sample with index i is given by the value at index i of each provided vector.");
            }

            if (design.ColumnCount > yy.Count)
            {
                throw new ArgumentException($"A regression of the requested order requires at least {design.ColumnCount} samples. Only {yy.Count} samples have been provided.");
            }

            //return design.QR().Solve(yy).ToArray();

            //QR<double> qr = design.QR();
            QR<double> qr = DenseQREx.Create(design, designCopy, QRMethod.Thin);
            var x2 = builder2.SameAs(yy, storage.ColumnCount);
            qr.Solve(yy, x2);
            return x2.ToArray();
        }
    }

}
